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Abstract 



1-^ ' A computational approach to describe the energy relaxation of a high-frequency vibrational mode 

Ph: 

in a fluctuating heterogeneous environment is outlined. Extending previous work [H. Fujisaki, 



Y. Zhang, and J.E. Straub, J. Chem. Phys. 124, 144910 (2006)], second-order time-dependent per- 
turbation theory is employed which includes the fluctuations of the parameters in the Hamiltonian 
within the vibrational adiabatic approximation. This means that the time-dependent vibrational 
frequencies along an MD trajectory are obtained via a partial geometry optimization of the solute 
with fixed solvent and a subsequent normal mode calculation. Adopting the amide I mode of 
N-methylacet amide in heavy water as a test problem, it is shown that the inclusion of dynamic 

> ; 

^^ . fluctuations may significantly change the vibrational energy relaxation. In particular, it is found 



that relaxation occurs in two phases, because for short times (^ 200 fs) the spectral density appears 



C^^ ' continuous due to the frequency-time uncertainty relation, while at longer times the discrete nature 



of the bath becomes apparent. Considering the excellent agreement between theory and experi- 
ment, it is speculated if this behavior can explain the experimentally obtained biphasic relaxation 



kS I the amide I mode of N-methylacetamide. 



I. INTRODUCTION 

Vibrational energy relaxation (VER) is ubiquitous in chemistry and physics. When a 
chemical reaction or conformational change occurs, some vibrational modes may be excited 
or de-excited, resulting in nonequilibrium phenomena of VER (for example, see Ref. [l| for 
the possible role of VER in enzyme reactions). It is usually assumed that VER is much 
faster than chemical reactions or conformational changes,^ but this is not always the case.-i^ 
Recent progress in time-resolved spectroscop5»^'^>^'^ii2i>iiii^'i^ as well as numerous theoretical 
formulation3i^ii^ii^ii^'i^ii^'22i^i2Ii2^iMi2^i2^i22i2^i2£i22i^i^i2^i^i2^ have been trying to unravel the 
molecular origin of VER. 

To define the problem, we consider the case of the photoinduced n = l—>-0 VER of a high- 
frequency (~ 1000—2000 cm~^) vibrational mode in a polyatomic molecule. We partition 
the total Hamiltonian H as 

H = Hs + Hb + Hsb, (1) 

where the system Hamiltonian Hs describes the high-frequency mode us, the bath Hamil- 
tonian Hb comprises all remaining vibrational modes Ua of the molecule, and Hsb accounts 
for the coupling between the system and bath. The dynamics of the total system is described 
by the Liouville equation 

ih^ = [H,pit)] (2) 

with p{t) being the density operator for the system and bath. Employing fully quantum- 
mechanical formulations, various quantum-classical approaches or quasiclassical methods, a 
number of groups have described VER by calculating the time evolution of the total system. 
While this strategy is formally exact, it is also quite time-consuming and therefore limited 
to small molecules or model systems. Thus, various reduced density-matrix formulations for 
the system part in Eq. ([2]) have been pursued in the fields of quantum optics, condensed 
matter physics, and physical chemistry.— i2Ii2^ 

Considering VER, it is often the case that the coupling V = Hsb — {Hsb)b (where {■..)b 
denotes the bath average) is small enough to be described by low-order time- dependent 
perturbation theory. Assuming, for simplicity, that the VER is dominated by cubic coupling 
Csa/3 between the system mode qs and two bath modes qa,(lp (i-e., V oc J^aB^Sa/sQsQalp), 



we obtain Fermi's Golden Rule^ii^ 

ki^o(x'^\{f\V\i)f5{uJs-uJa-(^p) (3) 

Q!,/3 

for the VER rate from the vibrationally excited state |1) to the ground state |0), where 
\i) and |/) denote the initial and final state of the total system, respectively. If we are 
only interested in the short-time decay of our initially prepared state (rather than the time 
evolution of the complete system as in Eq. ([2])), Fermi's Golden Rule provides a convenient 
and correct description of VER. However, in many cases (e.g., in condensed phase VER) the 
straightforward implementation of the Golden Rule is hampered by the fact that the true 
form of the bath Hamiltonian Hb and the coupling Hamiltonian Hsb is not known. As a 
remedy, either idealized models (e.g., a harmonic bath with bilinear couplings to the system 
mode) or classical molecular dynamics (MD) simulations have been invoked, which allow 
us to directly calculate the spectral density of the bath and its (inverse) Fourier transform, 
the bath correlation function {V{t)V{0)). Rewriting Eq. ([3]) in its time- dependent form, the 
Golden Rule can be directly expressed in terms of the bath correlation function 

/oo 
dte''^^f^'{V{t)V{0)), (4) 

-oo 

where Aufi accounts for the frequency difference of the initial and final state. In the case 
of high-frequency modes, a classical bath correlation function (sampled at ksT) may be 
only a poor approximation to the quantum correlation function (dominated by zero point 
energy motion). Hence it is a well-established approach to augment Eq. (jlj) with quantum 
correction factors.— 'i^'^^'i^ Another extension of Eq. (jlj) was given by Bakker, who explicitly 
included the time-dependent fluctuations of the bath frequencies^. 

Recently, Fujisaki et a/.^^'^^ proposed an alternative approach to include in a realistic 
manner the effects of an inhomogeneous environment into the description of VER. The idea 
is to first select numerous snapshots of the solute molecule and its surrounding solvent shell 
from an equilibrium MD trajectory. For each structure, instantaneous normal mode analysis 
of the solvated system is performed, which gives the instantaneous vibrational frequencies 
Us and Ua of this conformation as well as the instantaneous vibrational couplings. Subse- 
quently, a time-dependent perturbative calculation (which can reduce to the Golden Rule) 
is performed for each snapshot, and the overall rate is obtained by a direct average over 
all conformation-dependent rates. The approach rests on the assumptions (i) that instan- 



taneous normal mode calculations provide a correct representation of the time-dependent 
vibrational frequencies and (ii) that these frequencies are constant on the time scale of 
the VER process, thus rendering a direct inhomogeneous averaging appropriate. A similar 
treatment using the path-integral method has been proposed by Shiga, Okazaki and their 
coworkers.^'^ 

The goal of this work is to go beyond these assumptions by explicitly considering the time- 
dependence of instantaneous vibrational frequencies and couplings during the VER process. 
To this end, we include the time-dependent driving of the environment into the perturbative 



formulation of Refs. l32lJ33l . Adopting N-methylacetamide (NMA) in heavy water as a simple 
peptide model to study the VER of the amide I vibration, we discuss the effects of the 
various treatments of the fluctuations. Furthermore, we show that, within the adiabatic 
approximation underlying the approach, the time-dependent vibrational frequencies need to 
be obtained from geometry-optimized (rather than true instantaneous) molecular structures. 
The calculations nicely reproduce the subpicosecond VER of the amide I mode in NMA 
measured in recent transient infrared experiments.-'^'^ 



II. THEORY AND METHODS 

A. Perturbation theory: Previous formulation 

Before considering the time-dependent driving, it is helpful to first briefly summarize the 
derivation of the perturbative formulation of VER given by Fujisaki et al.—'^ The total 
Hamiltonian ([1]) is partitioned as H = Hq + V with 

Ho = Hs + Hb + {Hsb)b, (5) 

V = HsB- {Hsb)b, (6) 

where {■ ■ ■)b denotes the average over the bath degrees of freedom. We have defined the 
interaction Hamiltonian V such that {V)b = 0, which is appropriate for the perturbative 
treatment.— Employing time-dependent perturbation theory with respect to the coupling 
V, the system density operator ps{t) = TrBp(t) in second order can be written as 

P^it) = ^TrB IdhJJ dt.Usit) [Vih), [Vih),pm] Ulit). (7) 

Here V{t) = t/g(t)Vf/o(t) represent the coupling V in the interaction representation and we 
have introduced unperturbed propagators f/fe(t) = e"*^*^*/^ with k = 0, S, B. 

To be specific, in the following we assume that the system and the bath is described in 
the harmonic approximation 

Hs = ^f + ^«l. (8) 

H, = y('4 + 4A, (9) 



a 



where pi and qi {i = S, a) are momenta and positions of the ith vibrational mode with 
frequency Ui. Furthermore, we assume that the VER is dominated by cubic coupling between 
the system mode and two bath modes, resulting in 

HsB = -qsF = qs^ Csapqaqp, (10) 

which noticeably differs from a bilinear coupling commonly employed.— i^i^i^ We consider 
the case that at time t = the total density operator factorizes in system and bath operator, 
i.e., p(0) = P5(0)pb(0), and that initially the system is in its first excited state, i.e., ps{^) = 



|1)(1|. Combining Eqs. ([7]) - (fTOj) . the population probability of the amide I vibrational 
ground state can be written as 

pSw = {o\pf{tm 

= Al(0|g5|l)|' / dt, r dh [((5F(t2-ti)<5F(0))Be*"^(*^-*^) + c.c.] (11) 
II' Jo Jo 

with 5F{t) = F{t)-{F)B and F{t) = Ul{t)FUB{t). Assuming, for simphcity, that only bath 

modes with huja > ksT ~ 200 cm~^ are of importance in the VER of amide I vibrations, 

the bath correlation function can be written as22 

{6F{t2)5F{ti))B = — y ^^e-'(""+"^)(*^-*i) . (12) 

2 ^ uj^up 

Insertion into Eq. (ITT|) yields the final result for the time-dependent ground state population 

(^)u\- ^ s^ ^safs [1 - cos(u;g - ^g - u;^)t] . . 

To recover the corresponding Golden Rule expression ([3]), the long-time limit of dpoo/dt is 



taken." Following Ref. 



32l . however, here we define the VER probability as 



„(2)/ 



P{t) = e-fooit). (14) 

At short times t ^ 300 fs, P{t) ^ 1— Pqq (t). At longer times, Eq. (HM may be advantageous, 
since it avoids problems associated with the violation of positivity of the density matrix (see 
below) . 

Finally, we account for the inhomogeneity of the environment by evaluating the VER 
probability (fT4l) for a number of MD snapshots {r = 1, . . . , A^}. Each snapshot provides 
instantaneous vibrational frequencies tUg (t), Ua (t) and couplings Cg^p{t) of this confor- 
mation. The mean relaxation probability {P{t))s is then obtained by a statistical average 
over all conformation- dependent probabilities Pr{t) 

1 ^ 
{Pit))s = j;^Y.Pr{t). (15) 

Alternatively, we may approximate the mean relaxation probability by a second-order cu- 
mulant expansion, yielding 

{Pit))s = (l-[pS(t)+pS(t) + ...])^ 

^exp|-(pS(t)\ |. (16) 



6 



Figure [T] compares the outcome of the various possible definitions of the mean VER prob- 
abihty. [Results are obtained from partial energy minimization and dynamical averaging, see 
below.] We find that the direct averaging [Eqs. (TT^ and ([H])] and the cumulant approxima- 
tion [Eq. fITC]) ] give virtually the same results, while the naive ansatz {P(t))s = l — (pQQ{t) ) 
deviates after ^ 300 fs and eventually becomes negative. The latter is associated with the 
fact that the eigenvalues of a reduced density matrix can be negative.— Note that Figure 
[1] also motivates our ansatz (IT4l) for the VER probability of individual trajectories. As 
the cumulant approximation appears to provide the clearest derivation of the mean VER 
probability, in the following we will use this definition for {P(t))s- 




200 400 600 800 1000 

t(fs) 



FIG. 1: Mean vibrational energy relaxation probability {P{t))s as obtained for NMA in D2O. 
Compared is the outcome of direct averaging [Eqs. (fT5]) and ([H]), red line], cumulant approximation 
[Eq. (fT6]) . green line] and the naive ansatz {P{t))s = 1 — (Poo (*)) (blue line). 



B. New formulation: Time-dependent driving 



In the perturbation theory described above, the vibrational frequencies and coupling 
elements are assumed to be constant during the VER process.—"^ This is an assumption 
which may break down if the time scale of the VER is similar to that of the vibrational 
parameters. In this work, we extend this formulation by explicitly considering the time 
dependence of instantaneous vibrational frequencies us = ^s(t) and Ua = uja(t) and coupling 
elements Csaf3 = Csajsit) in the derivation of the population probability. As a consequence. 



the Hamiltonian of system and bath become time-dependent, Hs = Hs{t) and Hb = Hsit), 
thus yielding the corresponding propagators 

Us{t) = e-(^/^)/o^s(r)dr^ (^7^ 

Second-order time- dependent perturbation theory then leads to an expression that is for- 
mally quite similar to Eq. flTTl) . except that the time dependence of all vibrational parameters 
is considered. The resulting bath correlation function reads 

X exp <^ ~i / [LUair] + ujf3{T)]dT > , (19) 

and we obtain for the time- dependent ground-state population probability 

.(2)/4^ - ^ \^ f ^-^ f^ A-,. Csaf^{t2)Csa(3{tl) 



p)^^(t) = -> / dt2 I dti 

^t^-^O "^0 y^s{h)^sitl)^a{t2)UJp{ti)uJc,{t2)uJl3{ti) 

/•t2 

X COS / [usir] - ujair) - u;^(r)](ir . (20) 

Jtx 

The mean VER probability in Eq. flTSl) is again defined by averaging over all conformation- 
dependent probabilities obtained from various MD trajectories. Equation (!20l) represents 
the main theoretical result of this paper. We refer to this result as "dynamic treatment" 
of VER, because it takes into account the time-dependent fluctuations of the environment. 
Assuming constant frequencies and couplings, we of course recover Eq. flT^ . which can be 
regarded as the inhomogeneous limit of Eq. (I20l) . 

C. Calculation of time-dependent frequencies 

In the formulation of VER developed above, we need to calculate instantaneous vibra- 
tional frequencies oJiif) and couplings Csaisit) along an MD trajectory. As the VER of the 
amide I mode also involves the vibrations of the first water shell of the peptide, we include 
NMA and the 16 nearest water molecules in the normal mode analysis (see Sec. Ill Dl below: 
the effect of the size of the included solvent shell is discussed in Fig. E]). There are several 
ways to calculate the normal modes of this subsystem for an instantaneous molecular confor- 



mation. As adopted in Refs. 



32 



33l . for example, one may simply perform an instantaneous 



normal mode calculation at this molecular geometry. Alternatively, one may first perform a 
partial geometry optimization of NMA (i.e., keeping the 16 solvent water fixed), before the 
normal mode calculation is done. 






FIG. 2: MD snapshot of N-methylacetamide (NMA) and its 16 nearest water molecules. 

The latter procedure is sometimes referred to as "quenched normal modes" method,— 
and is based on the underlying adiabatic approximation for the calculation of the (high- 
frequency) amide I mode. To explain this, we consider the standard Born-Oppenheimer 
approach to molecular electronic structures. Here the nuclear (i.e., slow) degrees of freedom 
coordinates are fixed, while the Schrodinger equation for the electronic (i.e., fast) degrees of 
freedom is solved. In a direct analogy, to calculate the structure of high-frequency vibrational 
modes, we keep the slow solvent degrees of freedom fixed and solve the problem of the fast 
degrees of freedom. Within the harmonic approximation, the latter amounts to a geometry 
optimization of the subsystem followed by a normal mode calculation. 

D. Simulation details 

All simulations were performed using the CHARMM simulation program package. ^^ We 
employed the CHARMM22 all-atom force field^ to model the solute NMA (H3C-COND- 
CH3) and the TIP3P water model^^ with doubled hydrogen masses to model the solvent 
D2O. We also performed simulations for fully deuterated NMA (D3C-COND-CD3). The 
peptide was placed in a periodic cubic box of (25.5 A)^ containing 551 D2O molecules. All 
bonds containing hydrogen bonds were constrained by SHAKE algorithm^ with a relative 



geometric tolerance of 10~^. We used a 10 A cutoff with a switching function for the 
nonbonded interaction calculations. After a standard equilibration protocol, we ran a 100 
ps NVT trajectory at 300 K, from which 100 statistically independent configurations were 
sampled. For each initial condition, a 1 ps NVE run was performed for the VER calculations, 
using a time step of 1 fs. 

For the normal mode calculations, the Hessian matrix with respect to the mass-weighted 
Cartesian coordinates Xi 

^^^ = -7=ir^ (21) 

of the system NMA/(D20)i6 was calculated by using the vib command of CHARMM. We 
obtain in total 180 normal modes. In the partial minimization scheme, the cons fix sele 
command was employed to constrain all atoms of water except NMA. The min command 
was used to minimize the energy of the subsystem NMA. The cubic couplings Csaf3 with 
respect to the normal modes qi were calculated from the Hessian matrix using numerical 
differentiation^ 

C - 1 "^'^ 
■^"^ 2 dqsdqo^dqp 

^ -l^l^U^JJ,, ^^^- , (22) 

where {Uia} comprises the eigenvectors of the Hessian matrix. To avoid problems with 
low-frequency bath modes whose frequency may become zero (or even imaginary in the case 
of instantaneous normal modes) in the denominator of Eqs. flT3l) and fl20l) . we only include 
bath modes with uja > 100 cm~^. 
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III. COMPUTATIONAL RESULTS 



A. Time-dependent frequencies 



Let us first consider the time-dependent vibrational frequencies along a typical trajec- 
tory by performing instantaneous normal mode analysis of NMA and its 16 nearest water 
molecules. Figure [3] (top panels) shows the time evolution of the amide I mode (mode 141) 
and several other modes that are of importance in the VER of the amide I mode (see Table 
HI] below) . The vibrational frequencies obtained from the normal mode calculations at the 
instantaneous (not minimized) molecular structures (left panel) are seen to undergo sub- 
stantial fluctuations. The range of fluctuations appears quite unrealistic, when compared to 
typical experimental infrared line widths. The right panel of Fig. [3] shows results obtained 
from normal mode calculations for partially minimized molecular structures. We see that 
the partial geometry optimization of the subsystem NMA/(D20)i6 results in a reduction 
of the frequency fluctuations to a physically reasonable range. This finding indicates - in 
accordance with the underlying adiabatic approximation - that normal mode calculations 
should be performed after partial geometry optimization rather than at the instantaneous 
molecular structure. The same conclusion was reached in a recent study of the VER of 
isolated NMA on ab initio potential energy surfaces.— 

To assess the effect of the frequency fluctuations on the VER of the amide I mode, we next 
consider the time dependence of the resonance condition Ausa/sit) = ousit) — Uait) — ujjsit). 
As is evident from Eqs. flT^ and (120]) . the minima of \AuJsai3{t)\ largely determine the effi- 
ciency of the VER process. Adopting the most important mode combinations {S,a,(3} in 
the VER of the amide I mode of NMA (see Table [T] below) , Fig. [3] (middle panels) compares 
the time evolution of Aoosajsit) obtained from normal mode calculations with (right) and 
without (left) partial geometry optimization of the subsystem NMA. The considerable dif- 
ferences found for the two results of Ausajsit) again underlines the importance of using the 
appropriate (i.e., partially minimized) normal mode frequencies. 

Also shown in Fig. [3] are the cubic couplings Csapif) associated with the VER process. 
From their definition in Eq. fl22l) . it is clear that this quantity, too, depends on how the 
normal modes are calculated. As expected, we find that the results for Csa/sit) without 
minimization fiuctuate much more compared to the results after the partial minimization. 
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FIG. 3: Time evolution of the vibrational dynamics of NMA in D2O, as obtained from instantaneous 
normal mode analysis with (right) and without (left) partial energy minimization. Shown are 
(upper panels) the normal mode frequencies Wj of selected vibrational modes, (middle panels) the 
frequency mismatch Aujsa/sit) = u)s(t)—cOa{t)—oJi3{t) for several resonant bath mode combinations, 
and (lower panels) the corresponding third-order anharmonic couplings Csafsit)- 

We note that apart from the frequencies we also need the normal mode eigenvectors to 
calculate Csaisit)- For simplicity, however, we have assumed that the character of the modes 
does not change significantly during the VER process and therefore neglected the explicit 
time dependence of the eigenvectors. This assumption may break down, when it is applied 
to strongly fluctuating instantaneous normal modes. 
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B. Numerical evaluation of the dynamic VER formula 

We are now in a position to study the outcome of the various theoretical treatments of 
VER introduced above. To this end, Fig. H] shows the VER probabihty P{t) of the amide I 
mode for ten representative trajectories (thin hne) as well as the mean population averaged 
over 60 trajectories (thick line). We again compare results obtained from normal mode 
calculations with (right) and without (left) partial geometry optimization of the subsystem 
NMA. Furthermore, the dynamic treatment [Eq. (120|) ] as well as its inhomogeneous limit 
[Eq. flT3l) ] are considered. In all cases, we find three stages of the time evolution of P{t). 
First, we observe a quadratic initial decay for t < 10 fs, which can be directly derived as 
short-time expansion of Eq. flT^ . This is followed by an approximately linear decay for times 
up to ~ 200 fs. Corresponding to a frequency uncertainty of Auj ~ 50 cm~^, this is the 
time scale on which the normal mode spectrum {cuj} of the "total" system (NMA/(D20)i6) 
appears continuous, because Au '> ojs — uj^ — ojp. At longer times, the discrete nature of 
the bath becomes apparent, leading to dispersion of the population probability for various 
MD trajectories. 

As may be expected, the different time dependence of the normal mode frequencies with 
and without the partial minimization also results in different time evolutions of the VER 
probability. However, the stronger fluctuations of the frequencies and couplings obtained 
from the calculations without minimization (see Fig. [3]) do not necessarily result in en- 
hanced fluctuations of P{t). More important is whether the time-dependent fluctuations 
are taken into account in the perturbative calculation. After ^ 100 fs, we observe that the 
inhomogeneously approximated VER probabilities largely disperse, while the dynamic treat- 
ment leads to comparatively similar population decays for various MD trajectories. This 
is because in the former case we use a constant resonance condition ujs = ^a + ^/s, which 
may differ significantly for individual trajectories. In the dynamic treatment, on the other 
hand, the time integration in Eq. ( I20l) smoothens the time-dependent resonance condition 
ujsit) = uja(t) + ojpit). To roughly quantify the VER dynamics, we introduce an average 
time scale r via^ 



T ^2 ^ ^1 Jti d'^ ^2 ~ il 

with ti = 0.0 ps and ^2=0. 5 ps. (Quite similar results are obtained for ti = 0.5 ps and 
t2=l-0 ps.) As noted in the caption of Fig. |U we find that the relaxation time r as well as 
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FIG. 4: Amide I mode relaxation probability P{t) of NMA, as obtained from an instantaneous 
normal mode analysis with (right) and without (left) partial energy minimization. Shown are re- 
sults from (upper panels) the inhomogeneous averaging approximation and (lower panel) dynamics 
averaging. The VER times according to Eq. ([23]) are (a) r = 0.54 ± 0.23 ps, (b) r = 0.43 ± 0.15 
ps, (c) r = 0.93 ± 0.15 ps, and (d) r = 0.60 ± 0.16 ps. 

its fluctuation may depend significantly on the theoretical treatment of VER. 

It is instructive to compare the above results with calculations for fully deuterated N- 
methylacet amide. As shown in Fig. [5l deuteration significantly affects the relaxation. While 
there are again clear differences between the various theoretical treatments, the trends are 
not that clear as in the case of singly deuterated NMA shown in Fig. HI Since deuteration 
may change normal mode frequencies as well as vibrational couplings, this rises the question 
on the relaxation pathway underlying the VER of NMA. 

Following definition fl23l) of an average relaxation time scale r, the individual energy flow 
pathways can be characterized by introducing a path-specific VER decay times Tsafs- (Note 
that 1/r = J2ai3 l/''"5«/3-) Adopting again the sample trajectory shown in Fig. [3], Table U 
lists the dominant relaxation pathways characterized by the shortest tso/b- In the case of 
singly deuterated NMA, two mode combinations [(a,/5) = (108,116) and (109,116)] results 
in particularly fast relaxation. To further analyze the relaxation mechanism, we consider 



14 




200 400 600 800 1000 

t(fs) 





^V ' ' ''* 


0.8 


^v^ 


0.6 


>^^ 




^^^«^st- 


0.4 


^^x^sl 




N — ; 


0.2 


- 



200 400 600 800 1000 

t(fs) 




200 400 600 800 1OO0 

t(fs) 




200 400 600 800 1000 

t(fs) 



FIG. 5: Same as in Fig. H but for fully deuterated NMA. The VER times are (a) r = 0.89 ± 0.69 
ps, (b) T = 1.64 ± 0.32 ps, (c) r = 0.96 ± 0.28 ps, and (d) r = 1.09 ± 0.32 ps. 

TABLE I: Dominant energy flow pathways of singly deuterated NMA (upper panel) and fully 
deuterated NMA (lower panel). The vibrational energy relaxation from of the initially excited 
amide I mode (1^5) to the vibrational modes lo^ and bjp is characterized by their path-specific re- 
laxation times Tsafi (in ps), time-averaged Fermi resonance parameter {Fsa/3), frequency mismatch 
(AuJsaij) (in cm~-^), and 3rd order coupling element {Csap) (in kcal/mol/A), respectively. 



Mode a,/3 


109,116 


108,116 


109,115 


109,117 


110,116 


112,112 


108,117 


108,114 


109,114 


TSap 


3.3 


3.6 


6.5 


7.1 


17.5 


18.7 


21.0 


43.2 


44.7 


{FSc.(3) 


0.209 


0.176 


0.109 


0.043 


0.033 


0.046 


0.055 


0.030 


0.028 


{^^Saf}) 


9.5 


17.9 


14.6 


37.8 


164.6 


46.0 


21.5 


56.0 


37.0 


(Cscn) 


2.7 


4.2 


2.1 


2.3 


8.4 


3.1 


1.6 


2.2 


1.4 




Mode a,/3 


116,116 


109,116 


110,116 


108,116 


109,119 


081,109 


076,116 


079,116 


075,109 


TSaP 


24.6 


26.3 


33.5 


57.7 


78.7 


80.9 


107.2 


108.9 


142.8 


{Fscp) 


0.046 


0.054 


0.096 


0.017 


0.008 


0.002 


0.003 


0.002 


0.002 


(Auso^ii) 


252.3 


151.2 


60.1 


172.7 


94.6 


937.0 


565.4 


546.9 


974.7 


(Cscf,) 


18.9 


10.1 


7.7 


3.6 


0.93 


0.91 


1.2 


1.0 


0.9 



the time-averaged Fermi resonance parameter defined as^*^ 

{\Csap\) 



{FsalB) 




h 



LUSUJaUJfS / {\AuJSafi\)' 



(24) 
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where (...) denotes the time average over a single trajectory. We see that in most (but 
not all) cases a large Fermi resonance parameter (Fsap) also means a short corresponding 
relaxation time Tsa/s- Furthermore, it is interesting to note that it is the resonance factor 
AuJsafB rather than the anharmonic coefficient Csa/s that mostly determines -Fga/J and thus 
the relaxation pathway.— >^ 

In the case of fully deuterated NMA, on the other hand, there are no such strongly 
resonant modes. Instead we find that the subpicosecond decay of the overall VER probability 
P(t) is the result of numerous decay channels with relatively long path-specific decay times. 
The weak oscillations seen in Fig. Orefiect the detuning AusajS ~ 200 cm~^ of the two most 
important paths. Although quite similar in simulation and experiment ,-1^1^ interestingly, the 
mechanism of amide I VER is different for singly and fully deuterated NMA. 



TABLE II: Characterization of the normal modes that mainly participate in the vibrational energy 
relaxation of singly deuterated NMA (upper panel) and fully deuterated NMA (lower panel). Shown 
are the vibrational frequency coa and various projections on atomic coordinates, see Eq. (j25p . Mode 
141 is the amide I mode. 



Mode # 


107 


108 


109 


110 


111 


112 


113 


116 


117 


139 


140 


141 


(iLla) (cm^^) 


454 


571 


590 


752 


765 


864 


983 


1093 


1129 


1448 


1570 


1681 


Pco 


0.16 


0.25 


0.27 


0.39 


0.37 


0.00 


0.12 


0.15 


0.19 


0.36 


0.09 


0.92 


Pnd 


0.10 


0.51 


0.33 


0.28 


0.39 


0.63 


0.04 


0.11 


0.39 


0.16 


0.15 


0.04 


-PCH3{C) 


0.28 


0.18 


0.31 


0.27 


0.16 


0.00 


0.83 


0.32 


0.18 


0.22 


0.01 


0.04 


-PCH3(N) 


0.16 


0.05 


0.07 


0.05 


0.08 


0.36 


0.02 


0.42 


0.24 


0.26 


0.75 


0.00 


-Pwatcr 


0.31 


0.02 


0.02 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 




Mode # 


107 


108 


109 


110 


111 


112 


113 


116 


117 


139 


140 


141 


{uja) (cm^i) 


433 


549 


560 


651 


686 


787 


798 


963 


1010 


1293 


1485 


1675 


Pco 


0.08 


0.28 


0.29 


0.10 


0.17 


0.03 


0.19 


0.49 


0.07 


0.00 


0.46 


0.94 


Pnd 


0.04 


0.21 


0.40 


0.34 


0.20 


0.36 


0.23 


0.13 


0.12 


0.00 


0.39 


0.04 


PCH3(C) 


0.12 


0.40 


0.25 


0.41 


0.23 


0.37 


0.42 


0.30 


0.43 


0.00 


0.05 


0.02 


-PcHaCN) 


0.08 


0.10 


0.03 


0.14 


0.40 


0.23 


0.16 


0.07 


0.38 


0.00 


0.11 


0.00 


Avatcr 


0.69 


0.02 


0.03 


0.01 


0.00 


0.00 


0.00 


0.00 


0.00 


1.00 


0.00 


0.00 



To characterize the normal modes that mainly participate in the VER process, we have 
calculated their projection on various atoms of NMA and the surrounding water molecules. 
For example, the projection of normal mode a onto the CO bond is defined by 



P( 



CO 



E 



m 



(25) 



ieCObond 
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where the sum goes over the x, y, and z coordinates of the corresponding C and O atoms, 
respectively, and {Uia} comprises the eigenvectors of the Hessian matrix [see Eq. fl22|) ]. 
Listing these projections. Table [III shows that the resonant modes are mainly localized on 
NMA, rather than on the solvent water. We also notice that deuteration red-shifts some of 
the lower frequencies, as expected, but not the amide I mode frequency. This appears to be 
the main reason of the little resonance of the fully deuterated case. 
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FIG. 6: Influence of the solvent on the amide I mode relaxation in NMA. Including 16 (left), 32 
(middle), or all (right) water molecules in the calculations, the resulting overall vibrational energy 
relaxation appears to be quite similar. In the case of isolated NMA (thick line in left panel), 
however, the relaxation probability P{t) follows the solution-phase decay only up to ~ 300 fs. 
Note that in the absence of solvent molecules, all trajectories coincide. 

To further study the influence of the solvent. Fig. [6] compares the VER dynamics of NMA, 
when various numbers of D2O molecules are included in the normal mode and subsequent 
VER calculations. In all cases, we show the VER probability of 10 single trajectories for 
singly deuterated NMA, obtained from the dynamic VER treatment and partial energy 
minimization. Including 16 and 32 water molecules, the resulting amide I relaxation appears 
to be quite similar [Fig. [6] (a) and (b)]. As another test, we furthermore took into account the 
remaining 535 water molecules as flxed charges [Fig. [6] (c).] Although individual relaxation 
pathways may be different, the flgure indicates that the flrst solvation shell with 16 water 
molecules is sufficient to account for the flrst step of VER of the amide I mode. Also shown 
in the left panel of Fig. [6] is the case of isolated NMA. In the absence of solvent water, we 
flnd that the excited state population probability P{t) follows the solution-phase decay only 
up to ~ 300 fs, and deviates to a larger value at longer times. This is a consequence of 
the fact that while at short times a few resonant pathways dominate, there are numerous 
possibilities of "off-resonant pathways" that may contribute at longer times. 



17 



Finally we remark on the comparison of the above results with existing experimental data. 
For both singly and fully deuterated NMA, biphasic relaxation Pexp(^) = pe"*/"^^ + (1— p)e~*/'^2 
of the amide I mode has been reported.->^i^ The relaxation times are ti = 0.45 ps {p = 0.8) 
and r2 = 4 ps for the singly deuterated NMA,^ and ti = 0.38 ps {p = 0.5) and T2 = 2.1 
ps for the fully deuterated NMA.^ Comparing Pcxpit) with our results obtained from the 
dynamic VER treatment and partial energy minimization, Fig. [7] shows excellent agreement 
between theory and experiment - maybe better as it can be expected from a simple force 
field modeling of the normal modes of the system. Independent of force field uncertainties, 
however, is our finding that the VER of the amide I mode in NMA consists of two phases 
(see also Figs. [Hand E]). This is because for t ^ 200 fs the system's spectral density appears 
continuous due to the frequency-time uncertainty relation, while at longer times the discrete 
nature of the bath becomes apparent. Although it is tempting to speculate that this behavior 
can explain the experimentally obtained biphasic relaxation of NMA, it is at present not 
clear if the approximations involved in our description allow for this conclusion. (E.g., the 
validity of time- dependent perturbation theory deteriorates at longer times.) Nevertheless, 
the observation that amide I relaxation in larger peptides occurs in a monoexponential 
manner- is consistent with the fact that the above described effect is expected to vanish for 
larger peptides with higher spectral density. 
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t(fs) 
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FIG. 7: Comparison of calculated (red lines) and experimental (green lines) results for the VER 
of the amide I mode in NMA. Shown are data for singly (left) and fully (right) deuterated NMA. 
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IV. CONCLUDING REMARKS 

We have outlined a computational approach to describe the energy relaxation of a high- 
frequency vibrational mode in a fluctuating heterogeneous environment. Extending the pre- 
vious work,^^ we have employed second time-dependent perturbation theory, which includes 
the fluctuations of the parameters in the Hamiltonian within the vibrationally adiabatic 
approximation. This means that the time-dependent vibrational frequencies along an MD 
trajectory are obtained via a partial geometry optimization of the peptide with fixed solvent 
water and a subsequent normal mode calculation. Although it requires more computational 
effort, the partial geometry optimization is necessary, because its omission results in un- 
reahstically high fluctuations of the vibrational frequencies and other quantities (see Fig. 

ED. 

Adopting the amide I VER of NMA in heavy water as a test problem, we have shown 
that the inclusion of dynamic fluctuations may significantly change the time evolution of 
the VER probability P{t) (see Fig. Hj). After ^ 200 fs, we observe that the inhomogeneously 
approximated VER probabilities largely disperse, while the dynamic treatment leads to 
comparatively similar population decays for various MD trajectories. This is because in the 
former case we use a constant resonance condition us = uja + ^i3, while the time integration 
in the dynamic treatment averages over the time-dependent resonance condition. 

To characterize the dominant energy flow pathways of the amide I VER of NMA, we have 
introduced path-specific VER decay times and considered the Fermi resonance parameter of 
the relaxation process. In the case of singly deuterated NMA, mainly two combinations of 
bath modes were found to achieve the vibrational relaxation. In the case of fully deuterated 
NMA, on the other hand, we have not observed such strongly resonant modes. Instead we 
found that the subpicosecond decay of the overall VER probability P(t) is the result of 
numerous decay channels with relatively long path-specific decay times. In both cases, we 
observed that the VER of the amide I mode in NMA consists of two phases (see Fig. Hj). 
This is because for t ^ 200 fs the system's spectral density appears continuous due to the 
frequency-time uncertainty relation, while at longer times the discrete nature of the bath 
becomes apparent. Considering our excellent agreement between theory and experiment 
(see Fig. [7]), it may be speculated if this behavior can explain the experimentally obtained 
biphasic relaxation of NMA. 
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Two directions for further researches are apparent. On a short time scale (here ^ 200 
fs), our results suggest that a single normal mode calculation should be sufficient to de- 
scribe VER in peptides. The restriction to ultrashort time scales therefore allows us to 
employ high-quality ab initio methods to calculate the vibrational structure of even large 
molecules.—!^ On a longer time scale (here ^ 500 fs), on the other hand, the validity of 
time-dependent perturbation theory is expected to deteriorate. As a straightforward exten- 
sion, the vibrational dynamics may be described by stochastic Schrodinger equations^ or 
stochastic Liouville equations.— i^ 
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